Introduction

This tutorial provides a comprehensive introduction to a modeling and simulation workflow in Pumas. This tutorial will not get into the details of Pumas specifics, but instead provide a narrative on the lines of a regular workflow in our day to day work, with brevity where required to allow a broad overview. Wherever possible, cross-references will be provided to documentation and detailed examples that provide deeper insght into a particular topic.

As part of this workflow you will be introduced to various aspects such as

1. Data wrangling in Julia

2. Exploratory analysis in Julia

3. Continous and discrete data non-linear mixed effects modeling in Pumas

4. Model comparison routines, post-processing, validation etc.

The study and design

CTMNopain is a novel anti-inflammatory agent under preliminary investigation. A dose-ranging trial was conducted comparing placebo with 3 doses of CTMNopain (5mg, 20mg and 80 mg QD). The maximum tolerated dose is 160 mg per day. Plasma concentrations (mg/L) of the drug were measured at 0, 0.5, 1, 1.5, 2, 2.5, 3-8 hours.

Pain score (0=no pain, 1=mild, 2=moderate, 3=severe) were obtained at time points when plasma concentration was collected. A pain score of 2 or more is considered as no pain relief.

The subjects can request for remedication if pain relief is not achieved after 2 hours post dose. Some subjects had remedication before 2 hours if they were not able to bear the pain. The time to remedication and the remedication status is available for subjects.

The goal

We are expected to provide input for an optimal dose/dosing regimen to be carried forward in future trials. Two datasets are provided below, that can be downloaded using the links below.

1. pk_painscore.csv

2. pain_remedication.csv

using PumasTutorials
using Random
using CSV
using Pumas
using PlottingUtilities
using PumasPlots
using PumasApps
using PumasReports
using NCAUtilities
using Chain
using Latexify
using Dates
using PlotlyBase
using CairoMakie 
interactive!(false)

false

Data Wrangling

We start by reading in the two dataset and making some quick summaries.

Note: As a general convention during this example, I will refer to dataframes by ending the name of the variable with df and the Population version of that dataframe will be without the df to avoid confusion.

pkpain_df = DataFrame(CSV.File(joinpath(dirname(pathof(PumasTutorials)), "..", "data", "intro", "pk_painscore.csv"), missingstrings=["", "NA", "."]))
remed_df  = DataFrame(CSV.File(joinpath(dirname(pathof(PumasTutorials)), "..", "data", "intro", "pain_remedication.csv"), missingstrings=["", "NA", "."]))

Let's filter out the placebo data as we don't need that for the PK analysis.

pkpain_noplb_df = filter(x -> !(occursin.("Placebo", x.ARM)), pkpain_df)

1,440 rows × 6 columns

	ARM	ID	TIME	CONC	PAINRELIEF	DOSE
	String	Int64	Float64	Float64	Int64	Int64
1	A20_0_at2h	1	0.0	0.0	0	20
2	A20_0_at2h	1	0.5	1.15578	1	20
3	A20_0_at2h	1	1.0	1.37211	1	20
4	A20_0_at2h	1	1.5	1.30058	1	20
5	A20_0_at2h	1	2.0	1.19195	1	20
6	A20_0_at2h	1	2.5	1.13602	1	20
7	A20_0_at2h	1	3.0	0.873224	1	20
8	A20_0_at2h	1	4.0	0.739963	1	20
9	A20_0_at2h	1	5.0	0.600143	0	20
10	A20_0_at2h	1	6.0	0.425624	1	20
11	A20_0_at2h	1	7.0	0.363418	1	20
12	A20_0_at2h	1	8.0	0.304177	1	20
13	A80_0_at2h	2	0.0	0.0	0	80
14	A80_0_at2h	2	0.5	4.93492	1	80
15	A80_0_at2h	2	1.0	4.56849	1	80
16	A80_0_at2h	2	1.5	4.23436	1	80
17	A80_0_at2h	2	2.0	3.35444	1	80
18	A80_0_at2h	2	2.5	2.70046	1	80
19	A80_0_at2h	2	3.0	2.38586	1	80
20	A80_0_at2h	2	4.0	1.84095	1	80
21	A80_0_at2h	2	5.0	1.64535	1	80
22	A80_0_at2h	2	6.0	1.36486	1	80
23	A80_0_at2h	2	7.0	1.20802	1	80
24	A80_0_at2h	2	8.0	0.983679	1	80
⋮	⋮	⋮	⋮	⋮	⋮	⋮

do some data wrangling and plotting here

Non-compartmental analysis

Let's begin by peforming a quick NCA of the concentration time profiles and view the exposure changes across doses. The input data specicfication for NCA analysis () requires the presence of a route column and an amount column that specifies the dose. So, let's add that in.

#adding route variable
pkpain_noplb_df[!,:route] .= "ev"
# creating an `amt` column
pkpain_noplb_df[!,:amt] .= ifelse.(pkpain_noplb_df.TIME .== 0, pkpain_noplb_df.DOSE, missing)

Now, we map the data variables to the read_nca function that prepares the data for NCA analysis.

pkpain_nca = read_nca(pkpain_noplb_df,
                      id = :ID,
                      time = :TIME,
                      amt = :amt,
                      observations = :CONC,
                      group = [:DOSE],
                      route = :route)

A full NCA Report is now obtained for completeness purposes, but later we will only extract a couple of key metrics of interest.

pk_nca = run_nca(pkpain_nca, sigdig=3)
first(pk_nca.reportdf, 10)

10 rows × 40 columns (omitted printing of 31 columns)

	id	DOSE	doseamt	tlag	tmax	cmax	tlast	clast	clast_pred
	Int64	Int64	Int64	Float64	Float64	Float64	Float64	Float64	Float64
1	6	5	5	0.0	1.0	0.316115	8.0	0.110616	0.107079
2	8	5	5	0.0	0.5	0.335649	8.0	0.0962241	0.0960789
3	9	5	5	0.0	0.5	0.538566	8.0	0.0430812	0.0387015
4	12	5	5	0.0	1.0	0.215302	8.0	0.0664953	0.0653076
5	17	5	5	0.0	1.0	0.190482	8.0	0.0558857	0.0553108
6	24	5	5	0.0	1.0	0.264185	8.0	0.0246395	0.0236468
7	29	5	5	0.0	0.5	0.353126	8.0	0.0376847	0.0366889
8	32	5	5	0.0	0.5	0.311369	8.0	0.151186	0.143299
9	38	5	5	0.0	0.5	0.295493	8.0	0.0771599	0.0758813
10	40	5	5	0.0	0.5	0.44879	8.0	0.0650122	0.0645857

As CTMNopain's effect maybe mainly related to maximum concentration (cmax) or area under the curve (auc), we present some summary statistics using the summarize function from NCA.

strata = [:DOSE]
parms = [:cmax, :aucinf_obs]
output = summarize(pk_nca.reportdf; stratify_by = strata, parameters = parms)

6 rows × 9 columns (omitted printing of 2 columns)

	DOSE	parameters	extrema	geomean	geomeanCV	geostd	mean
	Int64	String	Tuple…	Float64	Float64	Float64	Float64
1	5	aucinf_obs	(0.9137, 3.39947)	1.53397	86.6964	1.32989	1.59819
2	5	cmax	(0.190482, 0.538566)	0.345132	374.621	1.29294	0.356087
3	20	aucinf_obs	(2.77495, 14.1139)	6.01975	23.4821	1.41356	6.3764
4	20	cmax	(0.933113, 2.70061)	1.43399	88.0787	1.26304	1.47354
5	80	aucinf_obs	(13.7102, 49.122)	28.2973	4.74071	1.34149	29.5023
6	80	cmax	(3.29685, 8.47195)	5.64125	22.2948	1.25771	5.78671

f = parameters_vs_group(pk_nca, :cmax)[1]

Dose normalized PK parameters, cmax and aucinf were essentially dose proportional between for 5 mg, 20 mg and 80 mg doses. Based on visual inspection of the concentration time profiles as seen below, CTMNopain exhibited monophasic decline, and perhaps a one compartment model best fits the PK data.

pkpain_noplb_plot_df = filter(x -> !(x.TIME .== 0), pkpain_noplb_df)
f = summary_observations_vs_time(pkpain_nca, columns=1, rows=3)[1]
#how to convert to log scale - cc Michael
f

Pharmacokinetic modeling

As seen from the plot above, the concentrations decline monoexponentially. We will evaluate both one and two compartment structural models to assess best fit. Further, different residual error models will also be tested.

We will use the results from NCA to provide us good initial estimates. The mean clearance is 3.29, the mean volume is 16.45 and a good initial estimate for absorption rate as obtained by $0.693/(tmax/4)$ is 3.85

Data preparation for modeling

PumasNDF requires the presence of evid and cmt columns in the dataset.

pkpain_noplb_df[!, :evid] .= ifelse.(pkpain_noplb_df.TIME .== 0, 1, 0)
pkpain_noplb_df[!, :cmt] .= ifelse.(pkpain_noplb_df.TIME .== 0, 1, 2)
pkpain_noplb_df[!, :cmt2] .= 1 # for zero order absorption

Further, observations at time of dosing, i.e., when evid = 1 have to be missing

pkpain_noplb_df[!, :CONC] .= ifelse.(pkpain_noplb_df.evid .== 1, missing, pkpain_noplb_df.CONC)

The dataframe is now converted to a Population using read_pumas. Note that both observations and covariates are required to be an array even if it is one element.

pkpain_noplb = read_pumas(pkpain_noplb_df,
                          id           = :ID,
                          time         = :TIME,
                          amt          = :amt,
                          observations = [:CONC],
                          covariates   = [:DOSE],
                          evid         = :evid,
                          cmt          = :cmt)

Now that the data is transformed to a Population of subjects, we can explore different models.

One-compartment model

pk_1cmp = @model begin
  @metadata begin
    desc = "One Compartment Model"
    timeu = u"hr"
  end
  @param begin
    "Clearance (L/hr)"
    tvcl ∈ RealDomain(lower = 0, init = 3.2)
    "Volume (L)"
    tvv  ∈ RealDomain(lower = 0, init = 16.4)
    "Absorption rate constant (h-1)"
    tvka ∈ RealDomain(lower = 0, init = 3.8)
    """
    - ΩCL
    - ΩVc
    - ΩKa
    """
    Ω    ∈ PDiagDomain(init = [0.04,0.04,0.04])
    "Proportional RUV"
    σ_p  ∈ RealDomain(lower = 0.0001, init = 0.2)
  end
  @random begin
    η ~ MvNormal(Ω)
  end
  @covariates begin
    "Dose (mg)" 
    DOSE
  end
  @pre begin
    CL = tvcl * exp(η[1])
    Vc = tvv * exp(η[2])
    Ka = tvka * exp(η[3])
  end
  @dynamics Depots1Central1
  @derived begin
    cp := @. Central/Vc
    """
    CTMx Concentration (ng/mL)
    """
    CONC ~ @. Normal(cp, abs(cp)*σ_p)
  end
end

PumasModel
  Parameters: tvcl, tvv, tvka, Ω, σ_p
  Random effects: η
  Covariates: DOSE
  Dynamical variables: Depot, Central
  Derived: CONC
  Observed: CONC

Before going to fit the model, let's evaluate some helpful steps.

1. Simulation to check appropriatness of data and model

simpk = simobs(pk_1cmp, 
              pkpain_noplb, 
              init_params(pk_1cmp))
f = sim_plot(pk_1cmp,
            simpk[1:8])[1]
f                          
# add overlay of observed data (cc Michael)

Our NCA based initial guess on the parameters seem to work well.

Lets change the initial estimate of a couple of the parameters to evaluate the senstitivty.

pkparam = (init_params(pk_1cmp)..., tvka=2, tvv = 10)

(tvcl = 3.2, tvv = 10, tvka = 2, Ω = [0.04 0.0 0.0; 0.0 0.04 0.0; 0.0 0.0 0
.04], σ_p = 0.2)

simpk = simobs(pk_1cmp,
               pkpain_noplb,
               pkparam)
f = sim_plot(pk_1cmp,
            simpk)[1]
f                  
#add observations cc Michael

Changing the tvka and decreasing the tvv seemed to make an impact and observations go through the simulated lines.

To get a quick ballpark estimate of your PK parameters, we can do a NaivePooled analysis. Below we test the NaivePooled approach

NaivePooled

pkfit_np = fit(pk_1cmp,
               pkpain_noplb,
               init_params(pk_1cmp),
               Pumas.NaivePooled(),
               omegas = (:Ω,))

coeftable(pkfit_np)

7 rows × 2 columns

	parameter	estimate
	String	Float64
1	tvcl	3.005
2	tvv	14.0867
3	tvka	45.6113
4	Ω₁,₁	NaN
5	Ω₂,₂	NaN
6	Ω₃,₃	NaN
7	σ_p	0.32999

The final estimates from the NaivePooled approach seem reasonably close to our initial guess from NCA, except for the tvka parameter. We will stick with our initial guess

One way to be cautious before going into a complete fiting routine is to evaluate the likelihood of the individual subjects given the initial parameter values and see if anyone pops out as unreasonable. There are a few ways of doing this:

• check the loglikelihood subject wise

• check if there any influential subjects

Below, we are basically checking if the initial estimates for any subject are way off that we are unable to compute the initial loglikelihood.

lls = []
for subj in pkpain_noplb
  push!(lls,loglikelihood(pk_1cmp,
                   subj,
                   pkparam,
                   Pumas.FOCE()))
end
hist(lls; bins = 10, normalization = :none, color = (:black, 0.5))

The distribution of the loglikelihood's suggest no extreme outliers.

Now that we have a good handle on our data, lets go ahead and fit a population model

pkfit_1cmp = fit(pk_1cmp,
                 pkpain_noplb,
                 pkparam,
                 Pumas.FOCEI(),
                 constantcoef = (tvka = 2,))

infer(pkfit_1cmp)

Asymptotic inference results using sandwich estimator

Successful minimization:                      true

Likelihood approximation:               Pumas.FOCE
Log-likelihood value:                     1231.188
Number of subjects:                            120
Number of parameters:         Fixed      Optimized
                                  1              6
Observation records:         Active        Missing
    CONC:                      1320              0
    Total:                     1320              0

-------------------------------------------------------------------
        Estimate           SE                     95.0% C.I.
-------------------------------------------------------------------
tvcl     3.1642          0.086619         [ 2.9944  ;  3.334   ]
tvv     13.288           0.27482          [12.749   ; 13.827   ]
tvka     2.0             NaN              [  NaN    ;   NaN      ]
Ω₁,₁     0.084948        0.011024         [ 0.063342;  0.10655 ]
Ω₂,₂     0.048569        0.0063502        [ 0.036122;  0.061015]
Ω₃,₃     5.5809          1.2194           [ 3.191   ;  7.9708  ]
σ_p      0.10093         0.0057198        [ 0.08972 ;  0.11214 ]
-------------------------------------------------------------------

Notice that tvka is fixed to 2 as we don't have a lot of information before tmax. From the results above, we see that the parameter precision for this model is reasonable.

Just to be sure, let's fit a 2-compartment model and evaluate

pk_2cmp = @model begin
  @param begin
    "Clearance (L/hr)"
    tvcl ∈ RealDomain(lower = 0, init = 3.2)
    "Central Volume (L)"
    tvv  ∈ RealDomain(lower = 0, init = 16.4)
    "Peripheral Volume (L)"
    tvvp ∈ RealDomain(lower = 0, init = 10)
    "Distributional Clearance (L/hr)"
    tvq  ∈ RealDomain(lower = 0, init = 2)
    "Absorption rate constant (h-1)"
    tvka ∈ RealDomain(lower = 0, init = 1.3)
    """
    - ΩCL
    - ΩVc
    - ΩKa
    - ΩVp
    - ΩQ
    """
    Ω    ∈ PDiagDomain(init = [0.04,0.04,0.04, 0.04, 0.04])
    "Proportional RUV"
    σ_p  ∈ RealDomain(lower = 0.0001, init = 0.2)
  end
  @random begin
    η ~ MvNormal(Ω)
  end
  @covariates begin
    "Dose (mg)"
    DOSE
  end
  @pre begin
    CL = tvcl * exp(η[1])
    Vc = tvv *  exp(η[2])
    Ka = tvka * exp(η[3])
    Vp = tvvp * exp(η[4])
    Q = tvq * exp(η[5])
  end
  @dynamics Depots1Central1Periph1
  @derived begin
    cp := @. Central/Vc
    """
    CTMx Concentration (ng/mL)
    """
    CONC ~ @. Normal(cp, cp*σ_p)
  end
end

PumasModel
  Parameters: tvcl, tvv, tvvp, tvq, tvka, Ω, σ_p
  Random effects: η
  Covariates: DOSE
  Dynamical variables: Depot, Central, Peripheral
  Derived: CONC
  Observed: CONC

pkfit_2cmp = fit(pk_2cmp,
                 pkpain_noplb,
                 init_params(pk_2cmp),
                 Pumas.FOCEI(),
                 constantcoef = (tvka = 2,))

infer(pkfit_2cmp)

Asymptotic inference results using sandwich estimator

Successful minimization:                      true

Likelihood approximation:               Pumas.FOCE
Log-likelihood value:                    1827.5359
Number of subjects:                            120
Number of parameters:         Fixed      Optimized
                                  1             10
Observation records:         Active        Missing
    CONC:                      1320              0
    Total:                     1320              0

-------------------------------------------------------------------
        Estimate           SE                     95.0% C.I.
-------------------------------------------------------------------
tvcl     2.8451          0.086279         [ 2.676   ;  3.0143  ]
tvv     10.988           0.26139          [10.476   ; 11.501   ]
tvvp     5.0546          0.42074          [ 4.2299  ;  5.8792  ]
tvq      1.5282          0.12132          [ 1.2904  ;  1.7659  ]
tvka     2.0             NaN              [  NaN    ;   NaN      ]
Ω₁,₁     0.10042         0.013818         [ 0.073339;  0.1275  ]
Ω₂,₂     0.061796        0.0081345        [ 0.045852;  0.077739]
Ω₃,₃     1.1334          0.18237          [ 0.77593 ;  1.4908  ]
Ω₄,₄     0.4641          0.13355          [ 0.20234 ;  0.72587 ]
Ω₅,₅     0.24145         0.045102         [ 0.15305 ;  0.32985 ]
σ_p      0.048443        0.0012226        [ 0.046046;  0.050839]
-------------------------------------------------------------------

The 2 compartment model has a much lower objective function compared to the 1 compartment. Lets compare the estimates from the 2 models.

@chain coeftable(pkfit_2cmp) begin
  leftjoin(coeftable(pkfit_1cmp), on = :parameter, makeunique = true)
  rename!(:estimate => :pk2cmp, :estimate_1 => :pk1cmp)
end

11 rows × 3 columns

	parameter	pk2cmp	pk1cmp
	String	Float64	Float64?
1	tvcl	2.84515	3.16418
2	tvv	10.9884	13.288
3	tvvp	5.05455	missing
4	tvq	1.52816	missing
5	tvka	2.0	2.0
6	Ω₁,₁	0.100421	0.0849481
7	Ω₂,₂	0.0617956	0.0485685
8	Ω₃,₃	1.13336	5.58088
9	Ω₄,₄	0.464105	missing
10	Ω₅,₅	0.241452	missing
11	σ_p	0.0484426	0.10093

We perform a likelihood ratio test to compare the two nested models. The test statistic and the P-value clearly indicate that a 2 compartment model is better.

lrtest(pkfit_1cmp, pkfit_2cmp)

Statistic:          1190.0
Degrees of freedom:      4
P-value:               0.0

We should also compare the other metrics and statistics, such ηshrinkage, ϵshrinkage, aic, bic

@chain PumasReports._model_metrics(pkfit_2cmp) begin
  leftjoin(PumasReports._model_metrics(pkfit_1cmp), on = :Metric, makeunique = true)
  rename!(:Value => :pk2cmp, :Value_1 => :pk1cmp)
end

11 rows × 3 columns

	Metric	pk2cmp	pk1cmp
	String	Float64	Float64?
1	Estimation Time	9.46	2.74
2	LogLikelihood ($LL$)	1830.0	1230.0
3	$-2LL$	-3660.0	-2460.0
4	AIC	-3640.0	-2450.0
5	BIC	-3580.0	-2420.0
6	(η-shrinkage) $η₁$	0.038	0.0158
7	(η-shrinkage) $η₂$	0.0469	0.0402
8	(η-shrinkage) $η₃$	0.511	0.733
9	(η-shrinkage) $η₄$	0.266	missing
10	(η-shrinkage) $η₅$	0.19	missing
11	(ϵ-shrinkage) $CONC$	0.185	0.105

We next generate some goodness of fit plots to compare which model is performing better. To do this, we first inspect the diagnostics of our model fit.

res_inspect_1cmp = inspect(pkfit_1cmp)
res_inspect_2cmp = inspect(pkfit_2cmp)

gof_1cmp = goodness_of_fit(res_inspect_1cmp)[1]
gof_1cmp

gof_2cmp = goodness_of_fit(res_inspect_2cmp)[1]
gof_2cmp

These plots clearly indicate that the 2 compartment model is a better fit compared to the one compartment model. We can look at selected sample of individaul plots.

#rand_subjs = rand(1:length(pkpain_noplb), 9)
f = subject_fits(res_inspect_2cmp,
             separate = true, columns = 3, rows = 3)

We look at the first set of 9 individuals here by indexing into generated plot

f[6]

There a lot of important plotting functions you can use for your standard model diagnostics. Please make sure to read the documentation for plotting and the tutorial associated with it. Below, we are checking the distribution of the empirical Bayes estimates.

f = empirical_bayes_dist(res_inspect_2cmp)[1]

f = empirical_bayes_vs_covariates(res_inspect_2cmp, categorical = [:DOSE])[1]
# increase the height of this plot with resolution scaling cc Michael

Clearly,our guess at tvka seems off-target. Let's try and estimate tvka instead of fixing it to 2

pkfit_2cmp_unfix_ka = fit(pk_2cmp,
                 pkpain_noplb,
                 init_params(pk_2cmp),
                 Pumas.FOCEI())

infer(pkfit_2cmp_unfix_ka)

Asymptotic inference results using sandwich estimator

Successful minimization:                      true

Likelihood approximation:               Pumas.FOCE
Log-likelihood value:                    1899.2177
Number of subjects:                            120
Number of parameters:         Fixed      Optimized
                                  0             11
Observation records:         Active        Missing
    CONC:                      1320              0
    Total:                     1320              0

-----------------------------------------------------------------
        Estimate           SE                     95.0% C.I.
-----------------------------------------------------------------
tvcl     2.6384          0.10616          [ 2.4303  ;  2.8465  ]
tvv     11.36            0.26879          [10.834   ; 11.887   ]
tvvp     8.1962          1.0842           [ 6.0713  ; 10.321   ]
tvq      1.3182          0.076977         [ 1.1673  ;  1.4691  ]
tvka     4.8575          0.33975          [ 4.1916  ;  5.5234  ]
Ω₁,₁     0.1292          0.019136         [ 0.091697;  0.16671 ]
Ω₂,₂     0.060381        0.0074235        [ 0.045832;  0.074931]
Ω₃,₃     0.40714         0.083975         [ 0.24255 ;  0.57173 ]
Ω₄,₄     0.14067         0.061444         [ 0.020237;  0.26109 ]
Ω₅,₅     0.25354         0.048529         [ 0.15843 ;  0.34866 ]
σ_p      0.048809        0.0011716        [ 0.046513;  0.051106]
-----------------------------------------------------------------

#compare_estimates(pkfit_2cmp,pkfit_2cmp_unfix_ka)

Let's revaluate the goodness of fits and η distribution plots.

Not much change in the general gof plots

res_inspect_2cmp_unfix_ka = inspect(pkfit_2cmp_unfix_ka)
goodness_of_fit(res_inspect_2cmp_unfix_ka)[1]

But you can see a huge improvement in the ηka, (η₃) distribution which is now centered around zero

f = empirical_bayes_vs_covariates(res_inspect_2cmp_unfix_ka, categorical = [:DOSE])[1]
# increase the height of this plot with resolution scaling cc Michael

Finally looking at some individual plots for the same subjects as earlier

f = subject_fits(res_inspect_2cmp_unfix_ka,
             separate = true, columns = 3, rows = 3)

f[6]

The randomly sampled individual fits don't seem good in some individuals, but we can evaluate this via a vpc to see how to go about.

We can now perform a vpc to check.

pk_vpc = vpc(pkfit_2cmp_unfix_ka, 200; observations = [:CONC],
             stratify_by = [:DOSE],
             ensemblealg=EnsembleThreads())

Data Quantiles

99×4 DataFrame
 Row │ DOSE    time     CONC       τ
     │ Int64?  Float64  Float64    Float64
─────┼─────────────────────────────────────
   1 │     20      0.5  0.984084       0.1
   2 │     20      1.0  0.884154       0.1
   3 │     20      1.5  0.776766       0.1
   4 │     20      2.0  0.680191       0.1
   5 │     20      2.5  0.578424       0.1
   6 │     20      3.0  0.506193       0.1
   7 │     20      4.0  0.366619       0.1
   8 │     20      5.0  0.25662        0.1
  ⋮  │   ⋮        ⋮         ⋮         ⋮
  93 │      5      2.5  0.291732       0.9
  94 │      5      3.0  0.263513       0.9
  95 │      5      4.0  0.215232       0.9
  96 │      5      5.0  0.174649       0.9
  97 │      5      6.0  0.143855       0.9
  98 │      5      7.0  0.11426        0.9
  99 │      5      8.0  0.0940141      0.9
                            84 rows omitted
Simulation Quantiles

99×6 DataFrame
 Row │ τ        time     DOSE    lower      middle    upper
     │ Float64  Float64  Int64?  Float64    Float64   Float64
─────┼─────────────────────────────────────────────────────────
   1 │     0.1      0.5      20  0.869588   0.98296   1.08136
   2 │     0.1      1.0      20  0.799605   0.896778  0.972083
   3 │     0.1      1.5      20  0.720801   0.801354  0.86428
   4 │     0.1      2.0      20  0.637842   0.708761  0.77209
   5 │     0.1      2.5      20  0.559801   0.62364   0.68024
   6 │     0.1      3.0      20  0.477872   0.541218  0.600978
   7 │     0.1      4.0      20  0.336549   0.402695  0.454168
   8 │     0.1      5.0      20  0.231998   0.288964  0.339008
  ⋮  │    ⋮        ⋮       ⋮         ⋮         ⋮         ⋮
  93 │     0.9      2.5       5  0.279559   0.299476  0.322243
  94 │     0.9      3.0       5  0.253261   0.273235  0.295307
  95 │     0.9      4.0       5  0.209771   0.226646  0.24829
  96 │     0.9      5.0       5  0.168981   0.185933  0.20671
  97 │     0.9      6.0       5  0.135216   0.151848  0.171871
  98 │     0.9      7.0       5  0.107127   0.12366   0.145592
  99 │     0.9      8.0       5  0.0848253  0.101184  0.123274
                                                84 rows omitted

f = vpc_plot(pk_2cmp,
        pk_vpc, rows=1, columsns=1)

3-element Vector{Figure}:
 Figure()
 Figure()
 Figure()

f[1]

f[2]

f[3]

The visual predictive check suggests that the model captures the data well across all dose levels.

---